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Abstract 

Fractional calculus allows one to generalize the linear, one-dimensional, 
diffusion equation by replacing either the first time derivative or the second 
space derivative by a derivative of fractional order. The fundamental 
solutions of these generalized diffusion equations are shown to provide 
probability density functions, evolving on time or variable in space, which 
are related to the peculiar class of stable distributions. This property is 
a noteworthy generalization of what happens for the standard diffusion 
equation and can be relevant in treating financial and economical problems 
where the stable probability distributions are known to play a key role. 

1 Introduction 

Non-Gaussian probability distributions are becoming more common as data 
models, especially in economics where large fluctuations are expected. In 
fact, probability distributions with heavy tails are often met in economics 
and finance, which suggests to enlarge the arsenal of possible stochastic 
models by non-Gaussian processes. This conviction started in the early 
sixties after the appearance of a series of papers by Mandelbrot and 
his associates, who point out the importance of non-Gaussian probability 
distributions, formerly introduced by Pareto and Levy, and related scaling 
properties, to analyse economical and financial variables, as reported in 
the recent book by Mandelbrot (1997). Some examples of such variables 
are common stock prices changes, changes in other speculative prices, and 
interest rate changes. In this respect many works by different authors have 
recently appeared, see e.g. the recent books by Bouchaud & Potter (1997), 
Mantegna & Stanley (1998) and the references therein quoted. 

It is well known that the fundamental solution (or Green function) of 
the Cauchy problem for the standard linear diffusion equation provides at 
any time the probability density function (pdf) in space of the Gauss (or 
normal) law. This law exhibits all moments finite thanks to its exponential 
decay at infinity. In particular, the space variance of the Green function 
is proportional to the first power of time, a noteworthy property that 
can be understood by means of an unbiased random walk model for the 
Brownian motion, see e.g. Feller (1957). Less known is the property for 
which the fundamental solution of the Signalling problem for the same 
diffusion equation, provides at any position a unilateral pdf in time, known 
as Levy law, using the terminology of Feller (1966-1973). Because of its 
algebraic decay at infinity as i~ 3//2 , this law has all moments of integer 
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order divergent, and consequently its expectation value and variance are 
infinite. 

Both the Gauss and Levy laws belong to the general class of stable probability 
distributions, which are characterized by an index a (0 < a < 2), called 
index of stability or characteristic exponent. In particular, the index of the 
Gauss law is 2 , whereas that of the Levy law is 1/2 . 

In this paper we consider two different generalizations of the diffusion 
equation by means of fractional calculus, which allows us to replace either the 
first time derivative or the second space derivative by a suitable fractional 
derivative. Correspondingly, the generalized equation will be referred to 
as the time-fractional diffusion equation or the symmetric, space-fractional 
diffusion equation. Here we show how the fundamental solutions of this 
equation for the Cauchy and Signalling problems provide probability density 
functions related to certain stable distributions, so providing a natural 
generalization of what occurs for the standard diffusion equation. 

The plan of the paper is as follows. First of all, for the sake of convenience 
and completeness, we provide the essential notions of Riemann-Liouville 
Fractional Calculus and Levy Stable Probability Distributions in Appendix 
A and B, respectively. 

In Section 2, we recall the basic results for the standard diffusion 
equation concerning the fundamental solutions of the Cauchy and Signalling 
problems. In particular we provide the derivation of these solutions by the 
Fourier and Laplace transforms and the interpretation in terms of Gauss 
and Levy stable pdf , respectively. 

In Section 3, we consider the time-fractional diffusion equation and we 
formulate for it the basic Cauchy and Signalling problems to be treated in the 
subsequent two sections. Here we adopt the Riemann-Liouville approach to 
Fractional Calculus, and the related definition for the Caputo time- fractional 
derivative of a causal function of time. 

In Section 4, we solve the Cauchy problem for the time-fractional diffusion 
equation by using the technique of Fourier transform and we derive the 
corresponding fundamental solution in terms of a special function of Wright 
type in the similarity variable. In this case the solution can be interpreted 
as a noteworthy symmetric pdf in space with all moments finite, evolving 
in time. In particular, its space variance turns out to be proportional to a 
power of time equal to the order of the time-fractional derivative. 

In Section 5, we derive the fundamental solution for the Signalling problem 
of the time-fractional diffusion equation by using the technique of Laplace 
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transform. In this case the solution, still expressed in terms of a special 
function of Wright type, can be interpreted as a unilateral stable pdf in 
time, depending on position, with index of stability given by half of the 
order of the time-fractional derivative. 

In Section 6, we consider the symmetric, space-fractional diffusion equation. 
Here we adopt the Riesz approach to Fractional Calculus, and the related 
definition for the symmetric space-fractional derivative of a function of a 
single space variable. Here we treat the Cauchy problem by technique 
of Fourier transform and we derive the series representation of the 
corresponding Green function. In this case the fundamental solution is 
interpreted in terms of a symmetric stable pdf in space, evolving in time, 
with index of stability given by the order of the space-fractional derivative. 
To approximate such evolution we propose a random walk model, discrete 
in space and time, which is based on the Griinwald-Letnikov approximation 
of the fractional derivative. 

Finally, Section 7 is devoted to conclusions and remarks on related work. 



2 The standard diffusion equation 

For the standard diffusion equation we mean the linear partial differential 
equation 

d d 2 

— u(x,t)=V-^u(x,t), u = u(x,t), (2.1) 

where V denotes a positive constant with the dimensions L 2 T _1 , x and t 
are the space-time variables, and u = u(x, t) is the field variable, which is 
assumed to be a causal function of time, i.e. vanishing for t < . 
The typical physical phenomenon related to such an equation is the heat 
conduction in a thin solid rod extended along x , so the field variable u is 
the temperature. 

In order to guarantee the existence and the uniqueness of the solution, 
we must equip (1.1) with suitable data on the boundary of the space-time 
domain. The basic boundary- value problems for diffusion are the so-called 
Cauchy and Signalling problems. In the Cauchy problem, which concerns 
the space-time domain — oo < x < +oo , t > 0, the data are assigned at 
t = + on the whole space axis (initial data). In the Signalling problem, 
which concerns the space-time domain x > , t > , the data are assigned 
both at t = + on the semi-infinite space axis x > (initial data) and at 
x = + on the semi- infinite time axis t > (boundary data); here, as mostly 
usual, the initial data are assumed to be vanishing. 
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Denoting by g(x) and h(t) two given, sufficiently well-behaved functions, the 
basic problems are thus formulated as following: 

a) Cauchy problem 

■u(x,0 + ) = g(x) , -oo<x<+oo; u(=Foo,i) = 0, t>0; (2.2a) 

b) Signalling problem 

u(x,0 + ) = 0, x>0; u(0 + ,t) = h(t) , u(+oo,t) = 0, t>0. (2.26) 

Hereafter, for both the problems, we derive the classical results which will be 
properly generalized for the fractional diffusion equation in the subsequent 
sections. 

Let us begin with the Cauchy problem. It is well known that this initial value 
problem can be easily solved making use of the Fourier transform and its 
fundamental solution can be interpreted as a Gaussian pdf in x. Adopting 
the notation g(x) -j- g(n) with k £ R and 

/+oo . 
e +lKX g(x)dx, 
-oo 

1 f+OO 

g(x) = T' 1 [$(«)] = ^J oo e ~ lKX $(«) dn , 

the transformed solution satisfies the ordinary differential equation of the 
first order 

^- + K 2 V^j u( K ,t) =0, u(k,0+) =g(n), (2.3) 
and consequently it turns out to be 

u(M) =g{n)e~ K2vt . (2.4) 

Then, introducing 

g d c (x,t)^g^K,t) = e- K2vt , (2.5) 

where the upper index d refers to (standard) diffusion, the required solution, 
obtained by inversion of (2.4), can be expressed in terms of the space 
convolution u(x, t) = G%{£,, t) g(x — £) d£ , where 

gi M = -^—r^^l^Vt) . (2 . 6) 
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Here Q d {x,t) represents the fundamental solution (or Green function) of 
the Cauchy problem, since it corresponds to g(x) = S(x) . It turns out 
to be a function in x, even and normalized, i.e. Q^(x,t) = Q d (\x\,t) and 
G d (x, t) dx = 1 . We also note the identity 

\x\g d c (\x\,t) = ^M d ((), (2.7) 

where ( = \x\/(\/Vt l l 2 ) is the well-known similarity variable and 

M d (() = -^e~^ 2 / 4 . (2.8) 

We note that M d (() satisfies the normalization condition j °° M d (() d( = 1 . 



The interpretation of the Green function (2.6) in probability theory is 
straightforward since we easily recognize 



G d (x,t) =p G (x;a) :-- 



1 



2na 



e-*7( 2 0, a 2 = 2Vt, 



(2.9) 



where pc(x;a) denotes the well-known Gauss or normal pdf spread out 
over all real x (the space variable), whose moment of the second order, the 
variance, is a 2 . The associated cumulative distribution function (cdf) is 
known to be 



V G {x; a) := J p G (x'; a) dx' = - 1 + erf ) 



(2.10) 



where eif(z) := (2/y / 7r) Jq exp (— u 2 ) du denotes the error function. 
Furthermore, the moments of even order of the Gauss pdf turn out to be 
f+™x 2n p G {x;o-)dx = (2n- l)\la 2n , so 



/+oo 
x 2n g d {x,t)dx = (2n- \)\\{2Vt) n , n=l,2,.... (2.11) 
-oo 



Let us now consider the Signalling problem. This initial-boundary value 
problem can be easily solved by making use of the Laplace transform. 
Adopting the notation h{t) h(s) with s G C and 

h(s)=£[h(t)]= f°°e- st h(t)dt, 
Jo 
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h(t) = C' 1 \h(t)] = — I e st h(s)ds, 
L J 2iri J Br 

where Br denotes the Bromwich path, the transformed solution of the 
diffusion equation satisfies the ordinary differential equation of the second 
order 

/ d 2 s \ 

\Wv) fi (*' s )= ' u(0 + ,s) = h(s), u(+oo,s) = 0. (2.12) 
and consequently it turns out to be 

u(x,s) =h(s)e-( x /^ sl/2 . (2.13) 

Then introducing 

gf(x,t) + g*(x,s)=e-( x /^ sl/2 , (2.14) 

the required solution, obtained by inversion of (2.13), can be expressed in 
terms of the time convolution, u(x, t) = Jjj Q$(x, r) h(t — r) dr , where 

gd { t) = _^ t -3/2 e -x 2 /(4Vt) (215) 

Here Qg(x,t) represents the fundamental solution (or Green function) of the 
Signalling problem, since it corresponds to h(t) = 5(t) . We note that 

g«(x,t)= PL s(t^):=- / ^e-»/W, t>0, n = ^, (2.16) 



where PLs(t', A*) denotes the one-sided Levy-Smirnov pdf spread out over all 
non negative t (the time variable). The associated cdf is, see e.g. Feller 
(1966-1971) and Priiss (1993), 

V L (t; (j,) := f\ L {t'- M ) dt' = erfc (J^ = erfc (^=) , (2-17) 

where erfc (z) := 1 — erf (z) denotes the complenatary error function. 

The Levy-Smirnov pdf has all moments of integer order infinite, since it 
decays at infinity as t _3//2 . However, we note that the absolute moments of 
real order v are finite only ifO<z^<l/2.In particular, for this pdf the mean 
is infinite, for which we can take the median as expectation value. From 
T :> Ls(trned'i I 1 ) = 1/2, it turns out that t me d ~ 2/x , since the complementary 
error function gets the value 1/2 as its argument is approximatively 1/2. 
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We note that in the common domain x > , t > the Green functions of 
the two basic problems satisfy the identity 

xQ d c {x,t)=tg d s {x,t), (2.18) 

that we refer to as the reciprocity relation between the two fundamental 
solutions of the diffusion equation. Furthermore, in view of (2.7) and (2.18) 
we recognize the role of the function of the similarity variable, M d (Q , 
in providing the two fundamental solutions; we shall refer to it as to the 
normalized auxiliary function of the diffusion equation for both the Cauchy 
and Signalling problems. 



3 The time-fractional diffusion equation 

By the time-fractional diffusion equation we mean the linear evolution 
equation obtained from the classical diffusion equation by replacing the first- 
order time derivative by a fractional derivative (in the Caputo sense) of order 
a with < a < 2. In our notation it reads 

d a u d 2 u 

— =V — , u = u(x,t), 0<a<2, (3.1) 

where T> denotes a positive constant with the dimensions I? T~ a . From 
Appendix A we recall the definition of the Caputo fractional derivative of 
order a > for a (sufficiently well-behaved) causal function f(t) , see (A. 9), 

D% f(t) ■= r( 1 ; [\t - r) m - a /M(r) dr , (3.2) 
Y(m — a) Jo 

where m = 1, 2, . . . , and < m — 1 < a < m . According to (3.2) we thus 
need to distinguish the cases < a < 1 and 1 < a < 2 . In the the latter case 
(3.1) may be seen as a sort of interpolation between the standard diffusion 
equation and the standard wave equation. Introducing 

t x-i 

$a( * ):= F(a)' a>0 ' (3 ' 3) 

where the suffix + is just denoting that the function is vanishing for t < , 
we easily recognize that the equation (3.1) assumes the explicit forms : 
if < a < 1 , 
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if 1< a < 2 , 




/ (t-T) 
JO 



1-a 




Extending the classical analysis for the standard diffusion equation (2.1) to 
the above integro-differential equations (3.4-5), the Cauchy and Signalling 
problems are thus formulated as in equations (2.2), i.e. 

a) Cauchy problem 

u(x, + ) = g(x) , — oo < x < +oo ; u(^foo, t) = , t>0; (3.6a) 

b) Signalling problem 

u(x,0 + ) = 0, x>0; u{0 + ,t) = h(t) , u(+oo,t) = 0, t >0 . (3.66) 

However, if 1 < a < 2 , the presence in (3.5) of the second order time 
derivative of the field variable requires to specify the initial value of the first 
order time derivative ut(x,0 + ) , since in this case two linearly independent 
solutions are to be determined. To ensure the continuous dependence of our 
solution on the parameter a also in the transition from a = 1" to a = 1 + , 
we agree to assume ut(x, + ) = . 

We recognize that our fractional diffusion equation (3.1), when subject to 
the conditions (3.6), is equivalent to the integro-differential equation 



where < a < 2 . Such integro-differential equation has been investigated 
by several authors, including Schneider & Wyss (1989), Fujita (1990), Priiss 
(1993) and Engler (1997). 

In view of our subsequent analysis we find it convenient to put 



In fact the analysis of the time-fractional diffusion equation turns out to 
be easier if we adopt as a key parameter the half of the order of the 
time-fractional derivative. In future we shall provide the symbol a with 
other relevant meanings, as the index of stability of a stable probability 
distribution or the order of the space derivative in the space-fractional 
diffusion equation. 




(3.7) 



a 



< v < 1. 



(3.8) 
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Henceforth, we agree to insert the parameter v in the field variable, i.e. 
u = u(x, t; v) . By denoting the Green functions of the Cauchy and Signalling 
problems by G c (x, t; v) and G s (x, t; v) , respectively, the solutions of the two 
basic problems are obtained by a space or time convolution, u(x, t; v) = 
f-™Gc(€,t;v)g(x-£)d£, u(x,t;u) = $QQ s (x,T;v)h{t-T)dT , respectively. 
It should be noted that G c (x,t;v) = G c (\x\, t; v) , since the Green function 
turns out to be an even function of x . 

In the following two sections we shall compute the two fundamental solutions 
with the same techniques (based on Fourier and Laplace transforms) used 
for the standard diffusion equation and we shall provide their interpretation 
in terms of probability distributions. Most of the presented results are based 
on the papers by Mainardi (1994), (1995), (1996), (1997) and by Mainardi 
& Tomirotti (1995), (1997). 

4 The Cauchy problem for the time-fractional 
diffusion equation 

For the fractional diffusion equation (3.1) subject to (3.6a) the application 
of the Fourier transform leads to the ordinary differential equation of order 
a = 2z/, 



Using the results of Appendix A, see (A. 22-30), the transformed solution is 



where E2 V {-) denotes the Mittag-Leffler function of order 2v , and conse- 
quently for the Green function we have 



Since the Green function is a real and even function of x, its (exponential) 
Fourier transform can be expressed in terms of the cosine Fourier transform 
and thus is related to its spatial Laplace transform as follows 




(4.1) 




(4.2) 




(4.3) 




(4.4) 



g c (s,t;u) 



s=+ik 



+ G c {s,t; v) 



s=—ik 
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Indeed, a split occurs also in (4.3) according to the duplication formula for 
the Mittag-Leffler function, see (A. 26), 

g c (k,t;u) = E 2l/ (-K 2 Vt 2 ») = 

[E l/ {+inVvt u ) + E v {-inVVt v )}/2. 1 ' ' 

When v ^ 1/2 the inversion of the Fourier transform in (4.5) cannot be 
obtained by using a standard table of Fourier transform pairs; however, for 
any v £ (0, 1) such inversion can be achieved by appealing to the Laplace 
transform pair (A. 37) with r = \x\ , and s = ±in . In fact, taking into 
account the scaling property of the Laplace transform, we obtain from (4.5) 
and (A.37) 

g c (\x\,t-,u) = —^-M(-^--, V ), (4.6) 



2\/Vt u \\fVt- 

where M((; u) is the special function of Wright type, defined by (A. 31-33), 
and 

the similarity variable. We note the identity 

\x\g c (\x\,t;») = ^M((;v), (4.8) 

which generalizes to the time-fractional diffusion equation the identity (2.7) 
of the standard diffusion equation. Since M (£; v) d(, = 1, see (A. 40), 
the function M(Q;v) is the normalized auxiliary function of the fractional 
diffusion equation. 

We note that for the time-fractional diffusion equation the fundamental 
solution of the Cauchy problem is still a bilateral symmetric pdf in x (with 
two branches, for x > and x < , obtained one from the other by 
reflection), but is no longer of Gaussian type if v ^ 1/2 . In fact, for large 
\x\ each branch exhibits an exponential decay in the "stretched" variable 
l^li/C 1 -^) as can fog derived from the asymptotic representation (A. 36) of the 
auxiliary function M(-;v) . In fact, by using (4.7-8) and (A. 36), we obtain 

Gc(x, t; v) ~ o,(i) | x |("-V2)/(i-^) exp [_ 6 ^)| x |V(i-)] j (4.9) 

as \x\ — > oo , where a*(t) and 6*(t) are certain positive functions of time. 

Furthermore, the exponential decay in x provided by (4.9) ensures that all 
the absolute moments of positive order of g c (x, t; v) are finite. In particular, 
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using (4.8) and (A. 39) it turns out that the moments (of even order) are 

£j 2n Gc(x, t; v) dx = (Vt 2 T , n = 0,l,2,... (4.10) 

The formula (4.10) provides a generalization of the corresponding formula 
(2.11) valid for the standard diffusion equation, v = 1/2. Furthermore, we 
recognize that the variance associated to the pdf is now proportional to Vt 2u , 
which for v ^ 1/2 implies a phenomenon of anomalous diffusion. According 
to a usual terminology in statistical mechanics, the anomalous diffusion is 
said to be slow if < v < 1/2 and fast if 1/2 < v < 1 . 

In Figure 1, as an example, we compare versus \x\ , at fixed t, the 
fundamental solutions of the Cauchy problem with different v [y = 
1/4 , 1/2 , 3/4). We consider the range < \x\ < 4 and assume V = t = 1 . 



0.5 




Figure 1: The Cauchy problem for the time- fractional diffusion equation. 
The fundamental solutions versus \x\ with a) v = 1/4, b) v = 1/2, c) 
v = 3/4. 

We note the different behaviour of the pdf in the cases of slow diffusion [y = 
1/4) and fast diffusion {y = 3/4) with respect to the Gaussian behaviour 
of the standard diffusion (y = 1/2). In the limiting cases v = and v = 1 
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we have 



ftfetiO)-^, OA**!)- **-™**'*™ ■ ("I) 



We also recognize from the appendix B that for 1/2 < v < 1 any branch 
of the fundamental solution is proportional to the corresponding positive 
branch of an extremal stable pdf with index of stability a = \jv , which 
exhibits an exponential decay at infinity. In fact, applying (B.29) with 
a = 1/v and y = ( = \x\/{y/Vt v ) , from (4.7-8) we obtain 

Qd\x\,t;u) = — L-.l qi/v [|x|/(V©f); -(2-1/1/)' = 

(4.12) 

^■Pi/„(|z|; +1,1,0), Kl/i/<2. 
We also note that the stable distribution in (4.12) satisfies the condition 

r+oo 

J Pi/v (x; +1, 1, 0) dx = v , 1<1/i/<2. (4.13) 

5 The Signalling problem for the time-fractional 
diffusion equation 

For the fractional diffusion equation (3.1) subject to (3.6b) the application 
of the Laplace transform leads to the ordinary differential equation of order 
2, 

/ d 2 s 2v \ 

[ ^ - "p") u{x,s;v), u(0 + ,s;u) = h(s) , u(+oo, s; v) = . (5.1) 

Thus the transformed solution reads 

u{x,s;v) = ~h( s ) e -( x /^ sU , (5.2) 
so for the Green function we have 

G s (x, t; u) - g s (x, s; v) = ^1^) s " . (5.3) 
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When v ^ 1/2 the inversion of this Laplace transform cannot be obtained by 
looking in a standard table of Laplace transform pairs. Also here we appeal 
to a Laplace transform pair related to the Wright-type function M(£; v). In 
fact, using (A. 40) with r = t , and taking into account the scaling property 
of the Laplace transform, we obtain 

g a (x, t\v) = v — ^ M ( ■ v \ . (5.4) 

Introducing the similarity variable Q = x/(VVt u ) , we recognize the identity 

tg s (x,t;u) = v(M((;v), (5.5) 

which is the counterpart for the Signalling problem of the identity (4.8) valid 
for the Cauchy problem. 

Comparing (5.5) with (4.8) we obtain the reciprocity relation between the 
two fundamental solutions of the time-fractional diffusion equation, in the 
common domain x > , t > , 

2vxQ c (x, t; v) = t Q s (x, t; v) . (5-6) 



The interpretation of Q s (x,t;v) as a one-sided stable pdf in time is 
straightforward: in this respect we need to apply (B.28), with index of 
stability a = v and variable y = (^ 1 / u = t (\/V jx) x l v , in (5.5). We obtain 



G s (x,t; v) 



l/v 



x 



l/i/ 



X 



; —v 



Mt; -1,1,0). (5-7) 



In Figure 2, as an example, we compare versus t , at fixed x , the fundamental 
solutions of the Signalling problem with different v (y = 1/4, 1/2,3/4). We 
consider the range < t < 3 and assume V = x = 1 . 



We note the different behaviour of the pdf in the cases of slow diffusion 
(y = 1/4) and fast diffusion {y = 3/4) with respect to the Levy pdf for the 
standard diffusion [y = 1/2). In the limiting cases v = , 1 , we have 

g s (x,t;0) =5(t), g s (x,t;l) = 5(t-x/Vv). (5.8) 
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Figure 2: The Signalling problem for the time-fractional diffusion equation. 
The fundamental solutions versus t with a) v = 1/4, b) v = 1/2, c) 
v = 3/4. 



6 The Cauchy problem for the symmetric space- 
fractional diffusion equation 

The symmetric space-fractional diffusion equation is obtained from the 
classical diffusion equation by replacing the second-order space derivative by 
a symmetric space-fractional derivative (explained below) of order a with 
< a < 2 . In our notation we write this equation as 



d u d a u 
~dt ~ d\x\ a ' 



u = u(x,t;a) , x £ R, t £ Rq , 0<a<2, (6.1) 



where D is a positive coefficient with the dimensions L a T _1 . The 
fundamental solution for the Cauchy problem, Q c (x,t;a) is the solution of 
(6.1), subject to the initial condition u(x,0 + ;a) = 5(x) . 

The symmetric space-fractional derivative of any order a > of a sufficiently 
well-behaved function <p{x) , x £ R , may be defined as the pseudo- 
differential operator characterized in its Fourier representation by 



d a 
d\x\ a 



4>(x) -. , x,fceR, a>0. 



(6.2) 
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According to a usual terminology, — \n\ a is referred to as the symbol of our 
pseudo-differential operator, the symmetric space- fractional derivative, of 
order a . Here, we have adopted the notation introduced by Zaslavski, see 
e.g. Saichev & Zaslavski (1997). 

In order to properly introduce this kind of fractional derivative we need 
to consider a peculiar approach to fractional calculus different from the 
Riemann-Liouville one, already treated in Appendix A. This approach is 
indeed based on the so-called Riesz potentials (or integrals), that we prefer 
to consider later. 

At first, let us see how things become highly transparent by using an 
heuristic argument, originally due to Feller (1952). The idea is to start 
from the positive definite differential operator 

^ : =-^-* 2 = N 2 , (6-3) 

whose symbol is \n\ 2 , and form positive powers of this operator as pseudo- 
differential operators by their action in the Fourier-image space, i.e. 



+(MT =i«r «><>. (6.4) 




Thus the operator — A a l 2 can be interpreted as the required fractional 
derivative, i.e. 

A a / 2 = --f-, a>0. (6.5) 

d\x\ a 1 ' 

We note that the operator just defined must not be confused with a power 
of the first order differential operator ^ for which the symbol is —in . 

After the above considerations it is straightforward to obtain the Fourier 
image of the Green function of the Cauchy problem for the space-fractional 
diffusion equation. In fact, applying the Fourier transform to the equation 

(6.1) , subject to the initial condition u(x, + ; a) = 5(x) , and accounting for 

(6.2) , we obtain 

G c (x,t;a) = G c (\x\,t;a) + G c (k,t;a) = e~ Vt \ K \ a , 0<a<2. (6.6) 

We easily recognize that the Fourier transform of the Green function 
corresponds to the canonic form of a symmetric stable distribution of index 
of stability a and scaling factor 7 = (Vt) 1 ^ , see (B.8). Therefore we have 

G c (x,t;a) = Pa (x; 0, 7 ,0), 7 = (Vt) V« . (6.7) 
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For a = 1 and a = 2 we easily obtain the explicit expressions of the 
corresponding Green functions since in these cases they correspond to the 
Cauchy and Gauss distributions, 

M ^ = hi*Tm?' < 6 - 8) 

see (B.5), and 

g c ( X ,t;2))=' e-* 2 /^), (6.9) 
zy 7T u t 

in agreement with (2.6). 
We easily recognize that 

" = (6 - 10 > 

is the similarity variable for the space-fractional diffusion equation, in terms 
of which we can express the Green function for any a € (0, 2] . Indeed, we 
recognize that 

G c (x,t;a) = p t y /a gafaO), (6.H) 

where q a (rj;0) denotes the symmetric stable distribution of order a with 
Feller-type characteristic function, see (B. 14-15). Now we can express the 
Green function using the Feller series expansions (B. 21-22) with 9 = 0. We 
obtain: 

for < a < 1 , 

1 ^r(raa + l) / avr\ / - n \n 

g«fa°) = -— E n , sm n T H ) ' ( 6 - 12a ) 

' ?1=1 v / 

for 1< a < 2 , 

m=0 ^ ' 

In the limiting case a = 1 the above series reduce to geometrical series and 
therefore are no longer convergent in all of C . In particular, they represent 
the expansions of the function 01(77; 0) = l/[7r(l + 17 2 )] , convergent for r\ > 1 
and < i] < 1 , respectively. 

We also note that for any a £ (0,2] the functions g Q (?7;0) exhibit at the 
origin the value g Q (0;0) = T(l/a) /{it a) , and at the queues, excluding the 
Gaussian case a = 2 , the algebraic asymptotic behaviour, as r\ — > 00 , 

Qa(v,0) ~ ^ r (° + !) sin ( a |) ?T (a+1) , 0<a<2. (6.13) 
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In Figure 3, as an example, we compare versus x , at fixed t , the fundamental 
solutions of the Cauchy problem with different a (a = 1/2,1, 3/2,2). We 
consider the range — 6 < x < +6 and assume V = t = 1 . 



0.4 



0.2 
0.1 h 






a) 







0.4 



0.3 - 



0.2 



0.1 - 




6 -4 -2 2 4 6 -6 -4 -2 2 4 6 



Figure 3: The Cauchy problem for the simmetric space-fractional diffusion 
equation. The fundamental solutions versus x : plate a) a = 1/2 
(continuous line), a = 1 (dashed line); plate b) a = 3/4 (continuous line), 
a = 2 (dashed line). 

Let us now express more properly our operator (6.4) (with symbol \n\ a ) 
as inverse of a suitable integral operator I a whose symbol is M~ a . This 
operator can be found in the approach by Marcel Riesz to Fractional 
Calculus, see e.g. Samko, Kilbas k Marichev (1987-1993) and Rubin (1996). 



We recall that for any a > , a / 1,3,5,... and for a sufficiently well- 
behaved function <p(x) , x G R, the Riesz integral or Riesz potential I a and 
its image in the Fourier domain read 



I a (f>{x) :-- 



J— OO \rv\ 



2T(a) cos(vra/2) 



(6.14) 



On its turn, the Riesz potential can be written in terms of two Weyl integrals 
I± according to 



I a 4>{x) 



1 



2 cos(vra/2) 



[I«cf>(x)+I«cf>(x)] 



(6.15) 
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where 

I%cp{x) :=-L f X (x- O a_1 0(O^, 

1 ,+oc ( 6 - 16 ) 

/«<K*):=— - / (C-a:)"- 1 ^^. 

Then, at least in a formal way, the space-fractional derivative (6.2) turns 
out to be defined as the opposite of the (left) inverse of the Riesz fractional 
integral, i.e. 

' ( x ) ■= = 1 — k+X^) + lZ a (p{x)] . (6.17) 



<i|x| Q 2 cos(7ra/2) 

Notice that (6.14) and (6.17) become meaningless when a is an integer odd 
number. However, for our range of interest < a < 2 , the particular case 
a = 1 can be singled out since the corresponding Green function is already 
known, see (6.8). Thus, excluding the case a = 1, our space-fractional 
diffusion equation (6.1) can be re-written, x € R, t € Rq , as 

— = -Vr a u, u = u(x,t;a), 0< a < 2 , a / 1 , (6.18) 
where the operator I~ a is defined by (6.16-17). 

Here, in order to evaluate the fundamental solution of the Cauchy problem, 
interpreted as a probability density, we propose a numerical approach, 
original as far as we know, based on a (symmetric) random walk model, 
discrete in space and time, see also Gorenflo & Mainardi (1998a), Gorenflo 
& Mainardi (1998b) and Gorenflo, De Fabritiis & Mainardi (1999). We shall 
see how things become highly transparent, in that we properly generalize 
the classical random-walk argument of the standard diffusion equation 
to our space-fractional diffusion equation (6.18). So doing we are in 
position to provide a numerical simulation of the related (symmetric) stable 
distributions in a way analogous to the standard one for the Gaussian law. 

The essential idea is to approximate the left inverse operators I± a by the 
Griinwald-Letnikov scheme, on which the reader can inform himself in the 
treatises on fractional calculus, see e.g. Oldham & Spanier (1974), Samko, 
Kilbas & Marichev (1987-1993), Miller & Ross (1993), or in the recent review 
article by Gorenflo (1997). If h denotes a "small" positive step-length, these 
approximating operators read 

1 00 f \ 

hl± a <t>(x) := ^ E^ 1 )" [ k j H*Tkh). (6.19) 
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Assume, for simplicity, V = 1 , and introduce grid points xj = j h with 
h > , j £ Z , and time instances t n = nr with r > , n G No . Let there 
be given probabilities pjk > of jumping from point xj at instant t n to 
point Xfc at instant t n+ i and define probabilities yj(t n ) of the walker being 
at point Xj at instant t n . Then, by 

oo oo oo 

Vk{tn+l)= J2 Pj,kUj(t n ), J2 Ph k= J2 Pj,* = 1 ' ( 6 - 20 ) 
j=-oo fc=— oo i=— oo 

with pjfi = pkj , a symmetric random walk (more precisely a symmetric 
random jump) model is described. With the approximation 



Vi 



(tn) ~ / 



(xj+h/2) 



'(xj-h/2) 

and introducing the "scaling parameter" 



u(x, t n ) dx ~ h u(xj, t n ) , 



(6.21) 



(i 



we have solved 



n+l) 



h a 2|cos(avr/2)| ' 

-Vj(tn) _ 



fJ a yj(t n ), 



(6.22) 



(6.23) 



for yj(t n+ i) . So we have proved to have a consistent (for h — ► 0) symmetric 
random walk approximation to (6.18) by taking 
for < a < 1, < // < 1/2, 



h I- a yj (t n ) =fi^- [ h I+ a yj (t n ) + h lZ a yj (t n ) 
I Pjj = 1 - 2/x , Pj ,i±fc = m I Q I , > l ; 

«; for 1< a < 2 , < n < l/(2a) , 

/i a r 

hl~ a yj(tn) =L i — [hl+ a Vj+l(tn) + hl- a yj-l{tn) 
Pj,j = 1-2/1 a, p j>j±1 = fj, [1 + (2)] , 



(6.24) 



(6.25) 



(fc+i) 



k > 2, 



We note that our random walk model is not only symmetric, but also 
homogeneous, the transition probabilities Pj,j±k not depending on the index 
3 ■ 
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In the special case a = 2 we recover from (6.25) the well-known three-point 
approximation of the heat equation, because Pjj±k = f° r k > 2 . This 
means that for approximation of common diffusion only jumps of one step 
to the right or one to the left or jumps of width zero occur, whereas for 
< a < 2 (a ^ 1) arbitrary large jumps occur with power-like decaying 
probability, as it turns out from the asymptotic analysis for the transition 
probabilities given in (6.24-25). In fact, as k — ► oo , one finds 

Pj,j+k ~ r (« + !) sin («0 , < a < 2 . (6.26) 

This result thus provides the discrete counterpart of the asymptotic 
behaviour of the long power-law tails of the symmetric stable distributions, 
as foreseen by (6.13) when < a < 2 . 



7 Conclusions 

We have treated two generalizations of the standard, one-dimensional, 
diffusion equation, namely, the time-fractional diffusion equation and the 
symmetric space-fractional diffusion equation. For these equations we have 
derived the fundamental solutions using the transform methods of Fourier 
and Laplace, and exhibited their connections to extremal and symmetric 
stable probability densities, evolving on time or variable in space. For the 
symmetric space-fractional diffusion equation we have presented a stationary 
(in time), homogeneous (in space) symmetric random walk model, discrete 
in space and time, the step-lengths of the spatial grid and the time lapses 
between transitions properly scaled. In the limit of infinitesimally fine 
discretization this model (based on the Griinwald-Letnikov approximation 
to fractional derivatives) is consistent with the continuous diffusion process, 
i.e. convergent if interpreted as a difference scheme in the sense of numerical 
analysis^. 

From the mathematical viewpoint the field of such "fractional" general- 
izations is fascinating as there several mathematical disciplines meet and 
come to a fruitful interplay: e.g. probability theory and stochastic processes, 

2 Further generalizations have been considered by us and our collaborators in other 
papers, in which we have given a derivation of discrete random walk models related to 
more general space-time fractional diffusion equations. For a comprehensive analysis, see 
Gorenflo et al. (2002). Readers interested to the fundamental solutions of these fractional 
diffusion equations are referred to the paper by Mainardi et al. (2001) where analytical 
expressions and numerical plots are found. 
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integro-differential equations, transform theory, special functions, numerical 
analysis. As one may take from our References, one can observe that since 
some decades there is an ever growing interest in using the concepts of 
fractional calculus among physicists and economists. Among economists we 
like to refer the reader to a collection of papers on the topic of "Fractional 
Differencing and Long Memory Processes", edited by Baillie & King (1996). 

Appendix A: The Riemann-Liouville Fractional 
Calculus 

Fractional calculus is the field of mathematical analysis which deals with the 
investigation and applications of integrals and derivatives of arbitrary order. 
The term fractional is a misnomer, but it is retained following the prevailing 
use. This appendix is mostly based on the recent review by Gorenflo & 
Mainardi (1997). For more details on the classical treatment of fractional 
calculus the reader is referred to Erdelyi (1954), Oldham & Spanier (1974), 
Samko et al. (1987-1993) and Miller & Ross (1993). 

According to the Riemann-Liouville approach to fractional calculus, the 
notion of fractional Integral of order a (a > 0) is a natural consequence 
of the well known formula (usually attributed to Cauchy), that reduces the 
calculation of the n— fold primitive of a function fit) to a single integral of 
convolution type. In our notation the Cauchy formula reads 

J n f(t):=f n (t) = -^—[\t-TT- 1 f(T)dT, t>0, neN, (Al) 
(n — Ij! Jo 

where N is the set of positive integers. From this definition we note that 
f n {t) vanishes at t = with its derivatives of order 1,2, ...,n — 1. For 
convention we require that f(t) and henceforth f n (t) be a causal function, 
i.e. identically vanishing for t < 0. In a natural way one is led to extend 
the above formula from positive integer values of the index to any positive 
real values by using the Gamma function. Indeed, noting that (n — 1)! = 
r(n) , and introducing the arbitrary positive real number a , one defines the 
Fractional Integral of order a > : 

J a /(*):= -L (\t-Ty- 1 /(r)dr, t >0 , a G R+ , (A2) 
1 [a) Jo 

where R + is the set of positive real numbers. For complementation we define 
J° := I (Identity operator), i.e. we mean J° f(t) = f(t) . Furthermore, by 
J Q /(0 + ) we mean the limit (if it exists) of J a f(t) for t —> + ; this limit 
may be infinite. 
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We note the semigroup property J a J lS = J a+lS , a, > 0, which implies 
the commutative property J@ J a = J a J 13 , and the effect of our operators J a 
on the power functions 

JQ * 7= r(7 + t+ a) f7+Q ' a "°' 7>_1 ' t>0 ' (A3) 

These properties are of course a natural generalization of those known when 
the order is a positive integer. 

Introducing the ^Laplace transform by the notation C {f(t)} : = 
f£°e~ st f(t) dt = f(s) , s G C, and using the sign to denote a Laplace 
transform pair, i.e. f(t) -j- f(s) , we note the following rule for the Laplace 
transform of the fractional integral, 

J a f(t)^^-, a>0, (A4) 

which is the generalization of the case with an n-fold repeated integral. 

After the notion of fractional integral, that of fractional derivative of order 
a (a > 0) becomes a natural requirement and one is attempted to substitute 
a with —a in the above formulas. However, this generalization needs some 
care in order to guarantee the convergence of the integrals and preserve the 
well known properties of the ordinary derivative of integer order. 
Denoting by D n with n G N , the operator of the derivative of order n , 
we first note that D n J n = I , J n D n / I , n G N , i.e. D n is left-inverse 
(and not right-inverse) to the corresponding integral operator J n . In fact 
we easily recognize from (A.l) that 

n—l ,k 

J-D n f(t) = f(t)-Y / f {k) (0 + )-, t>0. (A.5) 

k=0 

As a consequence we expect that D a is defined as left-inverse to J a . For 
this purpose, introducing the positive integer m such that m — 1 < a < m , 
one defines the Fractional Derivative of order a > : 

D a f(t) := D m J m ' a f(t) , m-Ka<m, toGN, (A6) 

namely 

r d m r i /■* f(r) 



D a f(t)-. 



dt m 

d m 

I dt™ 



dr 



T(m - a) Jo (t - r ) a+1 - m 
f(t) , a = m 



m — 1 < a < m, 

(A.&) 
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Defining for compiementation D° = J° = I , then we easily recognize that 
D a J a = I, a > , and 

D a V= + a>0, 7 >-l, t>0. (A.7) 

r(7 + 1 - a) 

Of course, these properties are a natural generalization of those known when 
the order is a positive integer. 

Note the remarkable fact that the fractional derivative D a f is not zero for 
the constant function f(t) = 1 if a N . In fact, (A.7) with 7 = teaches 
us that 

D a l = * - r , a>0, t>0. (A8) 
L (1 — a) 

This, of course, is = for a £ N, due to the poles of the gamma function in 
the points 0,-1, —2, .... We now observe that an alternative definition of 
fractional derivative, originally introduced by Caputo (1967) (1969) in the 
late sixties and adopted by Caputo and Mainardi (1971) in the framework 
of the theory of Linear Viscoelasticity, is 

D° f(t) := J m ~ a D m f(t) m-l<a<m, mGN, (A9) 

namely 

r 1 ^ f( m) (T) 

— r / — 7— dr , m — 1 < a < m, 

D * a f(t) = I T(m - a) Jo {t - r)^ 1 m ^ 9/) 

- — f(t) , a = m. 

( dt m 

This definition is of course more restrictive than (A. 6), in that requires 
the absolute integr ability of the derivative of order m. Whenever we use 
the operator D% we (tacitly) assume that this condition is met. We easily 
recognize that in general 

D a f(t) := D m J m - a f(t) / J m ' a D m f(t) := f(t) , (A10) 

unless the function f(t) along with its first m — 1 derivatives vanishes at 
t = + . In fact, assuming that the passage of the m-derivative under the 
integral is legitimate, one recognizes that, for m — 1 < a < m and t > , 

m-l f k~a 

D a f{t) = D? f{t) + T(h _ ~n ' ( AU ) 

fc=0 1 ^ K a ^ L > 
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and therefore, recalling the fractional derivative of the power functions (A. 7), 

/ m-l f k \ 

D a [f{t) - Y, j^ f {k) (0 + )j=D-f(t). (A12) 

The alternative definition (A. 9) for the fractional derivative thus incorpo- 
rates the initial values of the function and of its integer derivatives of lower 
order. The subtraction of the Taylor polynomial of degree m — 1 at t = + 
from f(t) means a sort of regularization of the fractional derivative. In 
particular, according to this definition, the relevant property for which the 
fractional derivative of a constant is still zero can be easily recognized, i.e. 

£>?1 = 0, a>0. (A13) 

We now explore the most relevant differences between the two fractional 
derivatives (A. 6) and (A. 9). We agree to denote (A. 9) as the Caputo 
fractional derivative to distinguish it from the standard Riemann-Liouville 
fractional derivative (A. 6). We observe, again by looking at (A. 7), that 

D a t a-1 = 0, a > 0, t > 0. 

From above we thus recognize the following statements about functions 
which for t > admit the same fractional derivative of order a , with 
m — 1 < a < m , to G N , 

m 

D a f(t) = D a g(t) ^ f{t) = git) + £ Cj t^ , (A14) 

j'=i 

m 

m fit) = D« git) ^ fit) = git) + Y cj t m ~i . (A15) 

j'=i 

In these formulas the coefficients Cj are arbitrary constants. 

For the two definitions we also note a difference with respect to the formal 
limit as a — > (m — 1) + ; from (A. 6) and (A. 9) we obtain respectively, 

D a fit) -^D m J fit) = D m - 1 fit) ; (A16) 

D" fit) -^JD m fit) = D m - 1 fit) - /( m " 1 )(0+) . (A17) 

We now consider the Laplace transform of the two fractional derivatives. 
For the standard fractional derivative D a the Laplace transform, assumed to 
exist, requires the knowledge of the (bounded) initial values of the fractional 
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integral J m a and of its integer derivatives of order k = 1,2, . . . ,m — 1. The 
corresponding rule reads, in our notation, 

m— 1 

D a f(t) Hr s° /(a) - £ Z> fe J (m " a) /(0 + ) s" 1 - 1 -^ , (A18) 

fc=0 

where m — 1 < a < m . 

The Caputo fractional derivative appears more suitable to be treated by 
the Laplace transform technique in that it requires the knowledge of the 
(bounded) initial values of the function and of its integer derivatives of 
order k = 1, 2, . . . , m — 1 , in analogy with the case when a = m . In fact, by 
using (A. 4) and noting that 

r d« f(t) = r j m ~ a D m f(t) = j m D m f{t) = f(t) - J2 / (fc) (o + ) tj , 

fc=o K - 

(A19) 

we easily prove the following rule for the Laplace transform, 



m— 1 

£>" fit) ~7" s a fis) - f {k) (0 + ) s a - l " k , to - 1< a < m . (A.20) 

fc=0 

Indeed, the result (A. 20), first stated by Caputo (1969) by using the 
Fubini-Tonelli theorem, appears as the most "natural" generalization of the 
corresponding result well known for a = m . 

Gorenflo and Mainardi (1997) have pointed out the major utility of the 
Caputo fractional derivative in the treatment of differential equations of 
fractional order for physical applications. In fact, in physical problems, 
the initial conditions are usually expressed in terms of a given number 
of bounded values assumed by the field variable and its derivatives of 
integer order, no matter if the governing evolution equation may be a 
generic integro-differential equation and therefore, in particular, a fractional 
differential equatiorH- 

We now analyze the most simple differential equations of fractional order, 
including those which, by means of fractional derivatives, generalize the well- 
known ordinary differential equations related to relaxation and oscillation 



3 We note that the Caputo fractional derivative was so named after the book by 
Podlubny (1999). It coincides with that introduced, independently and a few later, 
by Dzherbashyan and Nersesyan (1968) as a regularization of the Riemann-Liouville 
fractional derivative. Nowadays, some Authors refer to it as the Caputo-Dzherbashyan 
fractional derivative. The prominent role of this fractional derivative in treating initial 
value problems was recognized in interesting papers by Kochubei (1989), (1990). 
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phenomena. Generally speaking, we consider the following differential 
equation of fractional order a > , 



/ m-l ,fc \ 

£>?u(t) = £>" (u(t) - Ejfei u(fc) ( 0+ ) J =-«(*) + ?(*)> *>°> ( A21 ) 

where u = u(t) is the field variable and q(t) is a given function. Here m is 
a positive integer uniquely defined by m - 1 < a < ra , which provides the 
number of the prescribed initial values u^ k \0 + ) = Ck , k = 0, 1, 2, . . . , m — 1 . 
Implicit in the form of (A. 21) is our desire to obtain solutions u(t) for which 
the (t) are continuous. In particular, the cases of fractional relaxation 
and fractional oscillation are obtained for < a < 1 and 1 < a < 2 , 
respectively 

The application of the Laplace transform through the Caputo formula (A. 20) 
yields 

m-l a-k-l i 

fi W=E*i^+T + F+Ta(*)- (A22) 

Now, in order to obtain the Laplace inversion of (A.22), we need to recall 
the Mittag-Leffler function of order a > , E a {z) . This function, so named 
from the great Swedish mathematician who introduced it at the beginning 
of this century, is defined by the following series and integral representation, 
valid in the whole complex plane, 



E a (z) = Y— —r = — — : / —da, a>0. (A23) 

^ o r(an + l) 2TTlJ H a<T a -Z 



Here Ha denotes the Hankel path, i.e. a loop which starts and ends at — oo 
and encircles the circular disk |<r| < \z\^ a in the positive sense. It turns out 
that E a (z) is an entire function of order p = 1/a and type 1 . 
The Mittag-Leffler function provides a simple generalization of the expo- 
nential function, to which it reduces for a = 1 . Particular cases from which 
elementary functions are recovered, are 

E 2 (+z 2 ) = cosh z , E 2 (-z 2 ) = cos z , z £ C , (A. 24) 

and 

E lj2 (±z x l 2 ) = e z \l + erf (±z 1 / 2 )] = e z erfc {^z 1 ' 2 ) , z G C , (A25) 
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where erf (erfc) denotes the (complementary) error function, defined as 

1 r z 2 



erf 



(z) := -= ( Z e~ u2 du , erfc (z) := 1 - erf (z) , zeC. 
yvr Jo 



A noteworthy property of the Mittag-Leffler function is based on the 
following duplication formula 

E a (z) = \ [E a/2 {+z l/2 ) + E^-z 1 ' 2 )] . (A26) 

In (A. 25-26) we agree to denote by z 1 ! 2 the main branch of the complex 
root of z . 

The Mittag-Leffler function is connected to the Laplace integral through the 
equation 

e - " E a (u a z) du = — !— a > . (A27) 
o 1 — z 

The integral at the L.H.S. was evaluated by Mittag-Leffler who showed that 
the region of its convergence contains the unit circle and is bounded by the 
line Re z l / a = 1 . The above integral is fundamental in the evaluation of the 
Laplace transform of E a (— At") with a > and A G C . In fact, putting in 
(A. 27) u = st and u a z = —Xt a with t > and A G C , we get the Laplace 
transform pair 



L 



s a-l 



E a (-\t a ) + ^-^, ResylX^. (A28) 
Then, using (A. 28), we put for k = 0, 1, . . . , m — 1 , 

u fc (t) := J k e a (t) - , e Q (t) := E a (-t a ) , (A29) 

and, from inversion of the Laplace transforms in (A. 22), we find 

m— 1 ,.i 



In particular, the formula (A. 30) encompasses the solutions for a = 1,2, 
since ei(i) = exp(— t) , e2(t) = cos t. When a is not integer, namely for 
m — l<a<m,we note that m — 1 represents the integer part of a 
(usually denoted by [a]) and m the number of initial conditions necessary 
and sufficient to ensure the uniqueness of the solution u(t). Thus the m 
functions Uk(t) = J k e a (t) with k = 0, 1, . . . , m — 1 represent those particular 
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solutions of the homogeneous equation which satisfy the initial conditions 
= &kh , h, k = 0, 1, . . . , m — 1 , and therefore they represent the 
fundamental solutions of the fractional equation (A. 21), in analogy with the 
case a = m . Furthermore, the function u§(t) = —u' (t) = —e' a (t) represents 
the impulse-response solution. 

The Mittag-Leffter function of order less than one turns out to be related 
through the Laplace integral to another special function of Wright type, 
denoted by M(z, v) with < v < 1 , following the notation introduced 
by Mainardi (1994, 1995). Since this function turns out to be relevant in 
the general framework of fractional calculus with special regard to stable 
probability distributions, we are going to summarize its basing properties. 
For more details on this function, see Mainardi (1997), Appendix A. 

Let us first recall the more general Wright function W\^(z),z € C, with 
A > — 1 and fi > . This function, so named from the British mathematician 
who introduced it between 1933 and 1941, is defined by the following series 
and integral representation, valid in the whole complex plane, 

where Ha denotes the Hankel path. It is possible to prove that the Wright 
function is entire of order 1/(1 + A) , hence of exponential type if A > . The 
case A = is trivial since Wo tfM (z) = e z /r( / u) . The case A = — v , fx = 1 — v 
with < v < 1 provides the function M(z, v) of special interest for us. 
Specifically, we have 

M{z-v):=W- u ^ v {-z) = ±-W- vfi {-z), 0<u<l, (A32) 
and therefore from (A. 31-32) 

1 00 f— z \n-l 

M(z;u) = — — Y(vn) sin (unir) 

*tA n - i y- (A33) 

\ { „ „rr v da 

= — / &°~ za —, — , 0<I/<1. 
2lTl J H a u 

In the series representation we have used the reflection formula for the 
Gamma function, T(x) T(l — x) = it/ sin irx . Explicit expressions of M(z; v) 
in terms of simpler known functions are expected in particular cases when 
v is a rational number. Relevant cases are u = 1/2, 1/3 for which 

M(z; 1/2) = -L exp (- z 2 /A) , (A34) 
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M(z; 1/3) = 3 2/3 Ai (z/3 1/3 ) , (A35) 
where Ai denotes the Airy function. 

When the argument is real and positive, i.e. z = r > 0, the existence of 
the Laplace transform of M(r; v) is ensured by the asymptotic behaviour, 
as derived by Mainardi & Tomirotti (1995), as r — > +oo , 



M{r/v; v) ~ a(i/) r(" ~ X / 2 )/ ^ ~ ") exp 



-b(u) r 



V(l 



(A36) 



where a(v) = l/y/2n(l - v) , b(y) = (1 - v)jv . 

It is an instructive exercise to derive the Laplace transform by interchanging 
the Laplace integral with the Hankel integral in (A. 33) and recalling the 
integral representation (A. 23) of the Mittag-Leffler function. We obtain the 
Laplace transform pair 

M(r; v) -T- E v [—s) , 0<u<l. (A.37) 

For v = 1/2, (A.37) with (A. 25) and (A. 34) provides the result, see e.g. 
Doetsch (1974), 

M(r; 1/2) := — exp (- r 2 /4) -=- E l/2 {-s) := exp (s 2 ) erfc (s) . (A.38) 



7T 

It would be noted that, since M(r, u) is not of exponential order, 
transforming term- by-term the Taylor series of M(r; v) yields a series of 
negative powers of s , which represents the asymptotic expansion of E u (—s) 
as s — > oo in a certain sector around the real axis. 

We also note that (A.37) with (A. 23) allows us to compute the moments of 
any real order S > of M(r; v) in the positive real axis. We obtain 



/ 

Jo 



-DC 



r S M(r;u)dr = ^ S + 1) S > . (A.39) 
1 [vo + 1) 



When 5 is integer we note that the moments are provided by the derivatives 
of the Mittag-Leffler function in the origin, i.e. 



r 

Jo 



-oo fin y( — I — 1 A 

r n M(r; v) dr = lim (-1)" — E v {-s) = , (AAO) 



s- 



where n = 0, 1,2,... . The normalization condition J* °° M (r; u) dr = 
E u (0) = 1 is recovered for n = 0. The relation with the Mittag-Leffler 
function stated in (A. 40) can be extended to the moments of non integer 
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order if we replace the ordinary derivative, of order n, with the corresponding 
fractional derivative, of order <5 ^ n, in the Caputo sense. 

Another exercise on the function M concerns the inversion of the Laplace 
transform exp(— s") , either by the complex integral formula or by the formal 
series method. We obtain the Laplace transform pair 

^M(l/rV) - exp(-0 , 0<u<l. (A41) 

For v = 1/2 , (A. 41) with (A. 34) provides the known result, see e.g. Doetsch 
(1974), 

^ M ^ 1/2 ; 1/2) := ^37, exp [- l/(4r)] * exp . (A42) 

We recall that a rigorous proof of (A. 41) was formerly given by Pollard 
(1946), based on a formal result by Humbert (1945). The Laplace transform 
pair was also obtained by Mikusihski (1959) and, albeit unaware of the 
previous results, by Buchen h Mainardi (1975) in a formal way. 



Appendix B: The Stable Probability Distributions 

The stable distributions are a fascinating and fruitful area of research in 
probability theory; furthermore, nowadays, they provide valuable models in 
physics, astronomy, economics, and communication theory. 

The general class of stable distributions was introduced and given this name 
by the French mathematician Paul Levy in the early 1920's, see Levy (1924, 
1925). The inspiration for Levy was the desire to generalize the celebrated 
Central Limit Theorem, according to which any probability distribution 
with finite variance belongs to the domain of attraction of the Gaussian 
distribution. 

Formerly, the topic attracted only moderate attention from the leading 
experts, though there were also enthusiasts, of whom the Russian 
mathematician Alexander Yakovlevich Khintchine should be mentioned first 
of all. The concept of stable distributions took full shape in 1937 with the 
appearance of Levy's monograph, see Levy (1937-1954), soon followed by 
Khintchine's monograph, see Khintchine (1938). 

The theory and properties of stable distributions are discussed in some 
classical books on probability theory including Gnedenko & Kolmogorov 
(1949-1954), Lukacs (1960-1970), Feller (1966-1971), Breiman (1968-1992), 
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Chung (1968-1974) and Laha & Rohatgi (1979). Also treatises on fractals 
devote particular attention to stable distributions in view of their properties 
of scale invariance, see e.g. Mandelbrot (1982) and Takayasu (1990). Sets of 
tables and graphs have been provided by Mandelbrot & Zarnfaller (1959), 
Fama & Roll (1968), Bo'lshev & Al. (1968) and Holt & Crow (1973). 

Only recently, monographs devoted solely to stable distributions and related 
stochastic processes have been appeared, i.e. Zolotarev (1983-1986), Janicki 
& Weron (1994), Samorodnitsky & Taqqu (1994), Uchaikin & Zolotarev 
(1999). We now can cite the paper by Mainardi, Luchko & Pagnini (2001) 
where the reader can find (convergent and asymptotic) representations and 
plots of the symmetric and non-symmetric stable densities generated by 
fractional diffusion equations. 

Stable distributions have three exclusive properties, which can be briefly 
summarized stating that they 1) are invariant under addition, 2) possess 
their own domain of attraction, and 3) admit a canonic characteristic 
function. 

Let us now illustrate the above properties which, providing necessary and 
sufficient conditions, can be assumed as equivalent definitions for a stable 
distribution. We recall the basic results without proof. 

A random variable X is said to have a stable distribution P(x) = Prob{X < 
x} if for any n > 2 , there is a positive number c n and a real number d n such 
that 

X 1 + X 2 + ... + X n = c n X + d n , (B.l) 

where X\, Xi, ■ ■ ■ X n denote mutually independent random variables with 

common distribution P(x) with X . Here the notation = denotes equality 
in distribution, i.e. means that the random variables on both sides have the 
same probability distribution. 

When mutually independent random variables have a common distribution 
[shared with a given random variable X], we also refer to them as 
independent, identically distributed (i.i.d) random variables [independent 
copies of X]. In general, the sum of i.i.d. random variables becomes 
a random variable with a distribution of different form. However, for 
independent random variables with a common stable distribution, the sum 
obeys to a distribution of the same type, which differs from the original 
one only for a scaling (c n ) and possibly for a shift (d n ). When in (B.l) the 
d n = the distribution is called strictly stable. 

It is known, see Feller (1966-1971), that the norming constants in (B.l) are 
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of the form 

c n = n 1 / a with 0<a<2. (B.2) 

The parameter a is called the characteristic exponent or the index of stability 
of the stable distribution. 

We agree to use the notation X ~ P a (x) to denote that the random variable 
X has a stable probability distribution with characteristic exponent a . We 
simply refer to P(x) , p(x) := dP/dx (probability density function = pdf) 
and X as a-stable distribution, density, random variable, respectively. 

The definition (B.l) with the theorem (B.2) can be stated in an alternative 
version that needs only two i.i.d. random variables, see also Lukacs (1960- 
1970). A random variable X is said to have a stable distribution if for any 
positive numbers A and B, there is a positive number C and a real number 
D such that 

AX 1 + BX 2 = CX + D, (B.3) 

where X\ and X 2 are independent copies of X . Then there is a number 
a e (0, 2] such that the number C in (B.3) satisfies C a = A a + B a . 

For a strictly stable distribution (B.3) holds with D = . This implies that 
all linear combinations of i.i.d. random variables obeying to a strictly stable 
distribution is a random variable with the same type of distribution. 

A stable distribution is called symmetric if the random variable —X has the 
same distribution. Of course, a symmetric stable distribution is necessarily 
strictly stable. 

Noteworthy examples of stable distributions are provided by the Gaussian 
(or normal) law (with a = 2) and by the Cauchy-Lorentz law (a = 1). The 
corresponding pdf are known to be 

po^a^J^-^e-^-^W), xeK, (BA) 

V27T<7 



where a 2 denotes the variance and [i the mean, and 

1 7 

pc{x;j,8) := - — iER, (B.5) 

TT (X — 0) Z + 7^ 

where 7 denotes the semi-interquartile range and 5 the " shift" . 

Another (equivalent) definition states that stable distributions are the only 
distributions that can be obtained as limits of normalized sums of i.i.d. 
random variables. A random variable X is said to have a domain of 
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attraction^, e. if there is a sequence of i.i.d. random variables Y±,Y2,... 
and sequences of positive numbers {^ n } and real numbers {S n }, such that 

Yl + Y2 + ... Yn+sAx 

In 

The notation 4> denotes convergence in distribution. 

It is clear that the previous definition (B.l) yields (B.6), e.g. , by taking the 
YiS to be independent and distributed like X . The converse is easy to show, 
see Gnedenko & Kolmogorov (1949-1954). Therefore we can alternatively 
state that a random variable X is said to have a stable distribution if it has 
a domain of attraction. 

When X is Gaussian and the Y{S are i.i.d. with finite variance, then (B.6) 
is the statement of the ordinary Central Limit Theorem. The domain 
of attraction of X is said normal when 7 n = n 1//a ; in general, 7 n = 
n x l a h(n) where h(x) , x > , is a slow varying function at infinity, that 
is, lim h(ux)/h(x) = 1 for all u > , see Feller (1971). The function 

x— >oo 

h(x) = logx , for example, is slowly varying at infinity. 

Another definition specifies the canonic form that the characteristic function 
(cf) of a stable distribution of index a must have. Recalling that the cf is 
the Fourier transform of the pdf, we use the notation p a (n) := (exp (inX)) 
p a (x) . We first note that a stable distribution is also infinitely divisible, i.e. 
for every positive integer n its cf can be expressed as the nth power of 
some cf. In fact, using the characteristic function, the relation (B.l) is 
transformed into 

[p a ( K )} n =p a (c nK )e id ^. (B.l) 

The functional equation (B.7) can be solved completely and the solution is 
known to be 

p a (n; (3, 7, 5) = exp {idn — 7° + i (sign k) (5oj(\k\, a)]} , (B.8) 

where 

u>{Ma) = {^?*{ 2 \> * a + \ (B.9) 
{ — (2/7t) log , 11 a = 1 . 

Consequently a random variable X is said to have a stable distribution if 
there are four real parameters a, {3, 7, <5 with < a < 2 , — 1 < f3 < +1 , 
7 > 0, such that its characteristic function has the canonic form (B.8-9). 
Then we write p a (x; (3, 7, 5) -^-p a (K; (3, 7, 5) and X ~ P a (x; (3, 7, 5) , so partly 
following the notation of Holt Sz Crow (1973) and Samorodnitsky & Taqqu 
(1994). 
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We note in (B.8-9) that (3 appears with different signs for a ^ 1 and a = 1 . 
This minor point has been the source of great confusion in the literature, see 
Hall (1980) for a discussion. The presence of the logarithm for a = 1 is the 
source of many difficulties, so this case has often to be treated separately. 

The cf (B.8-9) turns out to be a useful tool for studying a-stable distri- 
butions and for providing an interpretation of the additional parameters, 
(3 (skewness parameter), 7 (scale parameter) and 5 (shift parameter), see 
Samorodnitsky & Taqqu (1994). When a = 2 the cf refers to the Gaussian 
distribution with variance a 2 = 2 j 2 and mean \i = 5 ; in this case the value 
of the skewness parameter (3 is not specified because tan it = , and one 
conventionally takes (3 = 0. 

One easily recognizes that a stable distribution is symmetric if and only if 
(3 = 5 = and is symmetric about 5 if and only if (3 = . Stable distributions 
with extremal values of the skewness parameter are called extremal. One 
can prove that all the extremal stable distributions with < a < 1 are 
one-sided, the support being Rg~ if (3 = — 1 , and Rq if (3 = +1 . 

For the stable distributions P a (x;(3, 7,5) we now consider the asymptotic 
behaviour of the tail probabilities, T + (\) := Prob{X > A} and T~(X) := 
Prob {X < — A} , as A — > 00 . For the Gaussian case a = 2 the result is well 
known, see e.g. Feller (1957), 

1 e-^/(47 2 ) 

° = 2: TW ~V^ A ' A ^°°- (B - 10) 

Because of the above exponential decay all the moments of the corresponding 
pdf turn out to be finite, which is an exclusive property of this stable 
distribution. For all the other stable distributions the singularity of the 
characteristic function in the origin is responsible for the algebraic decay of 
the tail probabilities as indicated below, see e.g. Samorodnitsky & Taqqu 
(1994), 

0<a<2: lim A Q T ± (A) = C a j a (l^p)/2, (5.11) 

where 

/ f°° \ _1 f ~r \ 7 TT ' if a 7^ 1 , 

C a = I / x~ a sinxdx\ = <^ T(2 - a) cos (avr/2) {B.12) 

[ 2/tt , ifa = l. 

We note that for extremal distributions ((3 = ±1) the above algebraic decay 
holds true only for one tail, the left one if (3 = +1 , the right one if (3 = — 1 . 
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The other tail is either identically zero if < a < 1 (the distribution is 
one-sided !), or exhibits an exponential decay if 1 < a < 2 . Because of the 
algebraic decay we recognize that 

0<a<2: f Pa (x; (3,^,5) dx = 0{\~ a ), (5.13) 

J\x\>\ 

so the absolute moments of a stable non-Gaussian pdf turn out to be finite 
if their order u is < v < a and infinite if v > a . We are now convinced 
that the Gaussian distribution is the unique stable distribution with finite 
variance. Furthermore, when a < 1 , the first absolute moment is 
infinite as well, so we need to use the median to characterize the expected 
value. 

There is however a fundamental property shared by all the stable 
distributions that we like to point out: for any a the stable pdf are unimodal 
and indeed bell-shaped, i.e. their n-th derivative has exactly n zeros, see 
Gawronski (1964). 

We now come back to the cf of a stable distribution, in order to provide for 
a 7^ 1 and 5 = a simpler canonic form which allow us to derive convergent 
and asymptotic power series for the corresponding pdf. We first note that 
the two parameters 7 and S in (B.8), being related to a scale transformation 
and a translation, are not so essential since they do not change the shape 
of distributions. If we take 7 = 1 and 5 = , we obtain the so-called 
standardized form of the stable distribution and X ~ P a (x; (3, 1, 0) is referred 
to as the a-stable standardized random variable. Furthermore, we can choose 
the scale parameter 7 in such a way to get from (B.8-9) the simplified canonic 
form used by Feller (1952, 1966-1971) and Takayasu (1990) for strictly stable 
distributions (5 = 0) with a / 1 , which reads in an ad hoc notation, 

9a(«;<?) := J^y Ky Pa(y,e)dy = e W {-|«| Q e ±i07r / 2 } , (5.14) 

where the symbol ± takes the sign of k . This canonic form, that we refer to 
as the Feller canonic form, is derived from (B.8-9) if in addition to a ^ 1 
and 5 = we require 

y* = C os (0 , tan (0 I) = /3tan (a . (5.15) 

Here 9 is the skewness parameter instead of (3 and its domain is restricted 
in the following region (depending on a) 

\e\<{°> * °< a<1 > (5.16) 

1 1 - \ 2 - a , if l<a<2. v ' 
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Thus, when we use the Feller canonic form for strictly stable distributions 
with index a / 1 and skewness 9 , we implicitly select the scale parameter 
7 (0 < 7 < 1), which is related to a , (3 and 9 by (B.15). Specifically, the 
random variable Y ~ Q a (y;9) turns out to be related to the standardized 
random variable X ~ P a (x; (3, 1, 0) by the following relations 



with 



Y = X/ 7 , p a (x;(3,l,0)= 1 q a (y = 1 x;9), (B.17) 



l/a 



7 = [COS (07T/2)] 

6 = (2/tt) arctan [/3tan (an/ 2)] , , R 18 \ 

_ tan (g7r/2) 1 j 

tan (an/ 2) ' 

We recognize that q a (y,9) = q a (—y,—9), so the symmetric stable 
distributions are obtained if and only if 6* = . We note that for the 
symmetric stable distributions we get the identity between the standardized 
and the Levy canonic forms, since in (B.18) (3 = 9 = implies 7 = 1. 
A particular but noteworthy case is provided by p 2 (x; 0,1,0) = (72(2/; 0) , 
corresponding to the Gaussian distribution with variance o 2 = 2 . 

The extremal stable distributions, corresponding to (3 = ±1 , are now 
obtained for 9 = ±a if < a < 1 , and for 9 = =f(2 - a) if 1 < a < 2 ; for 
them the scaling parameter turns out to be 7 = [cos (\a\ 7T/2)] 1 /" . It may be 
an instructive exercise to carry out the inversion of the Fourier transform 
when a = 1/2 and 9 = —1/2 . In this case we obtain the analytical expression 
for the corresponding extremal stable pdf, known as the (one-sided) Levy- 
Smirnov density, 

qi/2 (y;-l/2) = ^ / =y-^e- 1 /( 4 y), y>0. (B.19) 

The standardized form for this distribution can be easily obtained from 
(B.19) using (B. 17-18) with a = 1/2 and 9 = -1/2. We get 7 = 
[cos (-vr/4)] 2 = 1/2 , p = -1 , so 

p 1/2 (x; -1, 1,0) = i q 1/2 (x/2; -1/2) = -L x^' 2 e" 1 / ( 2x ) , (£.20) 

^ V Z7T 



where x > , in agreement with Holt & Crow (1973) [§2.13, p. 147]. 
Feller (1952) has obtained from (B.14) the following representations by 
convergent power series for the stable distributions valid for y > , with 
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< a < 1 (negative powers), 



00 

E< 

n=l 



T(na + 1) . 
) — ; sm 



1 < a < 2 (positive powers), 



<?« = — E (-y) n - LJ —. — - 

n! 



sm 



2a 



a 



(£.21) 



fB.22) 



The values for y < can be obtained from (B. 21-22) using the identity 
Qa{—y, = Qa(y, —0) , y > . As a consequence of the convergence in all of 
C of the series in (B. 2 1-22) we recognize that the restrictions of the functions 
yia(y,d) on the two real semi-axis turn out to be equal to certain entire 
functions of argument l/|y| a forO<a<l and argument |y| for 1 < a < 2 . 



It has be shown, see e.g. Bergstrom (1952), Chao Chung-Jeh (1953), that the 
two series in (B. 21-22) provide also the asymptotic (divergent) expansions to 
the stable pdf with the ranges of a interchanged from those of convergence. 



From (B. 21-22) a relation between stable pdf with index a and 1/a can be 
derived as noted in Feller (1966-1971). Assuming 1/2 < a < 1 and y > 0, 
we obtain 

-^ TT qi/ a (y~ a ;0) = q a (y,e*), e* = a(9 + i)-i. (s.23) 

A quick check shows that 8* falls within the prescribed range, \9*\ < a, 
provided that \9\ < 2 — 1/a . 

We now consider two particular cases of the Feller series (B. 21-22), of 
particular interest for us, which turn out to be related to the entire function 
of Wright type, M(z; v) with < v < 1 , reported in Appendix A. These 
cases correspond to the following extremal distributions 

:=q a (y;-a), y>0, 0< a < 1 , (5.24) 

*2(V) :=q a (y;a-2), y>0, Ka<2, (B.25) 
for which the Feller series (B. 21-22) reduce to 

MV) = - jri-lT^y-™ 1 - 1 r(na + 1) sin(nTra) , y>0, (£.26) 
7r ^— ; n! 

n=l 
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and 

#2(9) = l£ ( M)»-^»-'IMfi±i) 8ill (=), y>0 . (B.27) 

7T 77/. \ / 

71=1 

In fact, recalling the series representation of the general Wright function, 
W\,n( z ) with A > — 1 , \x > , see (A. 31), and the definition of the function 
M(z; u) with < v < 1 , see (A. 32-33), we recognize that 

<Z> 1 (y) = ±W^ (-y- a ) = -^TM(y- a ;a), y >0 , (5.28) 

and 

*2(j/) = -W_ 1/Q)0 (-y) = ^M(y;l/a), y>0. (5.29) 

We would like to remark that the above relations with the Wright functions 
have been noted also by Engler (1997). 

It is worth to point out that, whereas <3?i(y) totally represents the one- 
sided stable pdf q a (y; —a) , < a < 1 , with support in TLq , $2(2/) is the 
restriction on the positive axis of q a (y; a — 2) , 1 < a < 2 , whose support is 
all of R . Since the function M(z; v) turns out to be normalized in Rj}", see 
(A. 39-40), we also note 

roc poo 

/ * 1 (y)dy = l; / $ 2 (j/)dy = l/a. (5.30) 

JO JO 



Using the results (A. 41) and (A. 37) we can easily evaluate the Laplace 
transforms of &i(y) and <J?2(y) j respectively. We obtain 

£[*i(y)] = $i(s) =exp(-s Q ), 0<a<l, (5.31) 

«] = ^) = ^i/ a H, Ka<2, (5.32) 

where Ey a (-) denotes the Mittag-Leffler function of order l/a , see (A. 23). 

It is an instructive exercise to derive the asymptotic behaviours of 3>i(y) and 
<3? 2 (y) as y —> + and y — > +00 . By using the expressions (5.28 — 29) in terms 
of the function M and recalling the series and asymptotic representations of 
this function, see (A. 33) and (A. 36), we obtain 

q ^-(2-a)/[2(l-a)] e - Cl ^/d-)^ ag y ^ Q+ ; 

'!',(//)-<; Q (5.33) 
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J/<i [l + 0(y)}, as y •()' 



$ 2 ( y ) = J T(l- 1/a) - ■ ' ( B 34) 



y (2-a)/[ 2 (a-l)] e - C22/ «/(- ag y ^ +00) 




where ci , C2 are positive constants depending on a . We note that the 
exponential decay is found for $i(y) as y — > + but as y — » +00 for $2(2/) • 
Explicit expressions for stable pdf can be derived form those for the function 
M(z; v) when v = 1/2 and v = 1/3, given in Appendix A, see (A. 34- 
35). Of course the v = 1/2 expression can be used to recover the well- 
known (symmetric) Gaussian distribution q2 (y; 0) accounting for (B.29), and 
the (one-sided) Levy distribution (71/2(2/; — 1/2), see (B.19), accounting for 
(B.28). The v = 1/3 expression provides, accounting for (B.28), 

Mi* -1/3) = 3- X/ V 4/3 Ai [(3y)-^\ =^y- 3/2 K 1 

^ (£.35) 
where Ai denotes the Airy function and Ki/ 3 the modified Bessel function of 
the second kind of order 1/3 . The equivalence between the two expressions 
in (B.35) can be proved in view of the relation, see Abramowitz & Stegun 
(1965-1972) [(10.4.14)], 

AlW-i^d^). (R36) 

The case a = 1/3 has also been discussed by Zolotarev (1983-1986), who 
has quoted the corresponding expression of the pdf in terms of K ]/ 3 . 

A general representation of all stable distributions (thus including the 
extremal distributions above considered) in terms of special functions has 
been only recently achieved by Schneider (1986). In his remarkable 
(but almost ignored) article, Schneider has established that all the stable 
distributions can be characterized in terms of a general class of special 
functions, the so-called Fox H functions, so named after Charles Fox (1961). 
For details on Fox H functions, see e.g. the books Mathai & Saxena (1978), 
Srivastava & AI. (1982) and the most recent paper by Kilbas and Saigo 
(1999). These functions are expressed in terms of special integrals in the 
complex-plane, the Mellin-Barnes integral^. 

4 The names refer to the two authors, who in the first 1910's developed the theory of 
these integrals using them for a complete integration of the hypergeometric differential 
equation. However, as pointed out in the the Bateman Project Handbook on High 
Transcendental Functions, see Erdelyi (1953), these integrals were first used by S. Pincherle 
in 1888. For a revisited analysis of the pioneering work of Pincherle (1853-1936, Professor 
of Mathematics at the University of Bologna from 1880 to 1928) we refer to the paper by 
Mainardi and Pagnini (2003). 
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